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Results to Dnte 

We have now received two Landsat-4 TM tapes (7 bands) within our study area in 
the southern Sierra Nevada (path 41, row 35. 10 December 1982, and path 42. row 34. 
13 January 1983). During the reporting period we have worked on registration of the 
TM data to digital topographic data, on comparison of TM. MSS and NOAA meteorologi- 
cal satellite data for snowcover mapping, and on radiative transfer models for atmos- 
pheric correction. 

These investigations are described in more detail below. 


••Account ledgers for June are incomplete. 
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Geography, Urilmqi, China. 
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9. May 5, 1983, Elmissvuity of Snow, Xinjiang Institute of Geography. Urilmqi. China. 
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Detailed Descriptions of Investigations 


1. EFraCT OF SPATIAL R^LUTION ON THE ACCURACY OF MAPPING SNOW COVERED 
AREA 


1.1. Synnpais 

Satellite data from three sensors, Landsat-4 TM and MSS, and NOA.A-7 AVUKR, will 
be used to map snow covered area (SCA) over selected basins to determine the error 
associated with using coarser spatial resolution. This information is needed to quantify 
the information lost when high resolution data are not available and to determine the 
minimum resolution requirements for basins of different sizes The effect of a discon- 
tinuous snowcover within a basin will also be evaluated. This effort will improve our 
ability to make quantitative comparisons of surface measurements made with different 
satellite instruments 

Satellite remote sensing has become increasingly important to hydrologists 
because it provides information on the spatial distribution of parameters of hydrologic 
importance. In snow and ice studies, remote sensing has been used to improve the 
monitoring of existing conditions and has been incorporated into several runoff fore- 
casting and management systems. Range [1979] presents a detailed review of 
research in this area through 1978, and the SVioie and Ice, chapter of the 5th Pecora 
Symposium [Deutsch et al., 1981] has descriptions of several applications currently 
implemented or being tested. 
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The most common operational use of remote sensing in snow studies is to monitor 
the areal extent of the snowcover [Bamas and Bovdey, 1974; Foster and Rango, 1975: 
Jl/cOtnnis et al., 1975; Rango and Itten, 1976; Rango and Salomtmson, 1976]. These 
efforts led other investigators to conclude that satellite derived estimates of snow 
covered area [SCA] could be used to improve forecasting of snowmelt runoff [Hanna- 
ford, 1977; Rango et at., 1977; Rango, 1978; Rango et at.. 1979; Brown et al., I960]. 
Martinec [1975] and Rango and Martinec [1979] carried this a step further in their 
model by including satellite derived measurements of SCA as an index to melt. Satel- 
lites can also be viewed as space-borne radiometers. Techniques have been developed 
to use these data to measure or estimate snow surface characteristics that are neces- 
sary for calculation of the surface radiation budget [Dozier, 1960a, 1961, 1962, 1963; 
Dozier et al., 1900, 1981; Dozier and Frew, 1981; Dozier and Warren. 1982; Frampton 
and Marks, I960; Marks, 1962]. While these e^orts go beyond mapping SCA, any use of 
satellite data for snow hydrology begins with the mapping cf the snow cover distribu- 
tion over the surface. To determine the spatiai accuracy of any type of satellite data, 
we begin with the estimate of SCA. 

The efforts cited above have all been i;ifluenced by the availability, cost, and 
characteristics of data from satellite instruments that are currently operational. 
Operational applications have usually been limited to the measurement of SCA from 
uncorrected photographic products using manual photo interpretation techniques. 
Digital snow mapping techniques, while faster and less subjective, have been limited to 
specific experiments because of processing difficulties and the high cost of satellite 
data in digital form. For these reasons no effort has been made to assess the relative 
spatial accuracy of utilizing satellite data in different forms or from different instru- 
ments to estimate SCA. The purpose of this investigation is to conduct an experiment 
to determine the advantages and disadvantages of using data from different satellite 
instruments. From this we will evaluate the effect of basin size and shape, and discort- 
tinuous snowcover on the accuracy of our estimate of SCA for several different types of 
satellite data. 

1.2. Satellite Data 

Environmental satellite data are available in a variety of forms that are useful to 
snow hydrologists. For our investigations on the effect of resolution we vrill use data 
from the Landsat Thematic Mapper (TM) and Multispectral Spectral Scanner (MSS), and 
from the NOAA Advanced Very High Resolution Radiometer (AVHRR). These are mul- 
tichannel instruments that are frequently used to map SCA over basins from lOkm^ to 
10,000 fcm 2. Table 1 presents a comparison of these instruments. 
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Table 1 




Characteristics erf TH, HSS, and AVHRR Satellite Data 



(spectral, spatial, temporal, and economic) 


Landsat'4 TTiematic Mapper (TM) 


ujouelengths 

nominal 

overpass 

digital data 

band 

(SOX amol. . um ) 

vixel size 

interval 

cost 

1 

.452 

.518 

28,5 m 

16 days 

S2500 

2 

.529 

.610 




3 

.624 

.693 




4 

.776 

.905 




5 

1.568 

1.784 




7 

2.097 

2.347 




6 

10.422 

11.661 

114 m 



Landsai-4 MuUispectral Scanner (MSS) 

1 

.497 

.607 

57 m 

16 days 

S650 

2 

.603 

.697 




3 

.704 

.814 




4 

.809 





NOAA-7 Advanced Very High Resolution Radwmeter (AVHRR) 

1 

.56 

.72 

1.1 km 

12 hrs. 

$75 

2 

.71 

.98 




3 

3.53 

3.94 




4 

10.32 

11.36 




l5 

11.45 

12.42 




Digital data costs are approximate for the TM as they are not yet available 
to the public. For this investigation we have the TM data in hand. Digital 
data costs for the AVHRR reflect our costs to travel to Scripps SSOF, use of 

the facility, and return to UCSB 

It is assumed that at least four data sets 

are acquired per trip. 





We have acquired TM and MSS data for the southern Sierra Nevada from NASA God- 
dard for December 10, 19B2. We have received an additional TM data set for the central 
Sierra Nevada for January 13. 19B3. We vrill acquire AVHRR data from the Scripps Satel- 
lite Oceanographic Facility (SSOF) for the same overpass dates, and we hope to acquire 
one additional MSS data set to correspond to the January TM overpass. 
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1.3. The Eiperiment 

By acquiring data from three different environmental satellites over the same area 
on the same date, we have the uncommon opportunity to compare the spatial accuracy 
associated with using each type to estimate SCA. We plan to answer the following fun* 
damental questions: 

« For basins of several different sizes and shapes, how much difference is there in 
our estimate of SCA when using TM. MSS, and AVHRR data? 

3R Under what conditions of basin size, and snow cover distribution will coarse resolu- 
tion AVHRR data provide adc.f,.,..le information on SCA? 

* How cem occasional inputs from TM or MSS data be used to improve estimates of 
SCA that are regularly made with AVHRR data? 

These questions are important because while high resolution TM data are undoubtedly 
more accurate, they are expensive, difficult to acquire with reasonable regularity, and 
beyond the digital processing capability of most researchers. MSS data are the most 
commonly used satellite data, but they also have poor temporal resolution, and are 
expensive. AVHRR data are inexpensive, and the satellite has frequent overpasses, but 
the coarse spatial resolution makes the data difficult to use in small basins, or under 
conditions of sparse or patchy snowcovers. 

The research is divided logically into two parts. The first step is to calibrate and 
register the three satellite data sets to a standard reference. This is by far the most 
difficult part of the work, and will take the most time and effort. The development of 
the QDIPS image processing system at the Computer Systems Laboratory. University of 
California, Santa Barbara, makes this registration possible. The second involves the 
selection of test basins of varying sizes emd snowcover distributions for the mapping of 
SCA using each type of satellite data. Each of these is described in detail below. 

1.3.1. Calibrating and Registering Satellite Data Sets fo a Reference Orid 
Before the satellite data can be analyzed, the following processing is required: 

» Radiometric correction: Pixel values must be converted to the appropriate physi- 
cal units (surface exitance or brightness temperature). 

* Geometric correction: The satellite magery must be geometrically transformed so 
that it will overlay our digital terrain database. 

The flow of processing is as follows: 

1.3.2. Radiometric correction 

TM, MSS. and AVHRR shortwave data are converted to surface exitances using a 
simple lineeu* treuisform whose multiplicative and additive terms are assumed to be 
constant [Lauritson et at., 1979; Engel, 1980; dark and Dasgupta, 1903]. Since these 
coefficients may change over an image frame time, the radiometric correction must be 
performed before any geometric corrections alter the original time-space relationships 
within the image. 

1.3.3. Geometric correction 

Our terrain database is constructed from USGS digital elevation data. The data 
are supplied in geodetic coordinates (latitude, longitude). During the process of 
extracting a region of interest, we transform the data to the Universal Transverse Mer- 
cator (UTM) projection. We have adopted the UTM projection as our geographic refer- 
ence since it maintains almost constant scale over areas of the sizes in which we are 
interested, eind is readily indexed in common units of measurement (meters). Accu- 
rate. computationally efficient formulas are available to convert between UTM and geo- 
detic (latitude and longitude) coordinates [Dozier, 1980b]. The spatial resolution of the 
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uulividiml It^iTAin ilcxlrt setn K i;rui!c*’) will be res(\mpleii lo ^ melers, 

eorrespoiuiu\|; lo the hi^ihesl satellile resolution O'W) to be used in this investigation 

Ovir problem is therefore to transform the TM. MSS. and A\ UKK d*ata into the ITM 
ooordinate system of the terrain ijrids. at the seleoted spatial resolution Koi l*M *and 
V.SS data, the standard digital piiHiuets supplied bv NASA in the Spav'e C'blu|ue V.erea 
tor ^Sl'V.^ projeotUMi Hie St'M pivjeetion lias been nuathematieallv defined SViyiier, 
but eomputer eode implementmi; the oonversions between St'M. and ^eodetie 
eoordinates is not eurrently available He will therefore have to develop programs for 
mapping from Si'V. into I’l'M Seale ehaiV|*es will be aeeommodated bv resampling Anv 
niisre|;istration will be eorreeted by polynomial fitting of ground eontrol points ^l.h'^s), 
b'xperienee of other users of TV. and MSS data leads us to believe that any sueh residual 
niisrei;istrations are lar^elv translational in nature, and tlierefore easilv eorreeted bv *a 
few thTs [/\iftd /Hff.v. NASA Johnson Spaee Kliijht I'enter. personal eommunieationl 

AVHKK data as obtained from the Satellite Oeeano^raphie Kaeilitv at the Seripps 
Institute of Oeeano^raphv have had no ^eometrie eorreetions performed Hie j^eneral 
problem of ^<eoreferenein^ satellite imaijerv requires preeise determin%'\tion of tlie 
position ^altitude and nadir point) and attitude ^vaw pit eh and roll) of the satellite, 
together with a model of the spatial resolutioi\ of the sensor svsteni He are eurrentlv 
implenientiivi; svstenis to prediet satellite positions vising ii fHvvftfrtarn v'rbital element 
data and to build a satellite attitude veetor from telemetered vaw piteh. and roll data 
lliese models are tuned by least -squares adjustment of their eoeftieients vising eontrol 
points with known ^eodetie loeations I'his piwedure will Ih* used to determine tlie 
mapping between AVHKK ima^e eoordinates and W\\ eoordinates. usii\|i eontrol points 
seleeted from pseudo linage displays of the terrain ^rid 

rhe end produet of the ^eometrie eorreetion proeess will be l‘M V.SS. and AVHK'K 
data sets registered to digital terrain nunlels and resampled lo kS meter resolution 
Sinee this resampling represents an interpolation of up to ‘^\\\ Hie AVHKK ease), we 
will investigate whether hi>;h order resampling teehniques afTeet SHA determination in 
any systematie fashion ywe suspeet that the efTeet will be ne^li^ible) 

I 4 Uap5»4My SVii>u\*oi erifd Art: a 

Several basins that are partiallv snoweovered and varv in si.-e from to 
will be seleeted To ftnd partiallv snoweowred basins that are small we will seleet areas 
near Hu^ Hwens vallev wl\ere small basins have a lai\<tn* elevation ehaiv^e Kstimates of 
Sc'A foi’ eaeh basin will be made vising eaeli satc'llite data type Hasins will be' di<i;itallv 
isolated usin^ a proeedure developed bv t/nrAw ef iii , ^ and mappiii^; of St'A will be 

done witl\ existing software in the ^HMKS system In eaeh ease' the TV data will be used 
as the stiifiAurvi to wietermine the best estimate of SC'A as well as the spa 

tial eonti^<uitv of the snoweover Kesults from the V.SS and A\HKl\ data will then be 
eompared to this stai\dard 

In addition to the diijital mappii\|; of St'A for eaeh test basin James Foster, an 
e\periei\eed photo ii\terpreter from NASA Goddard Spaee Fli>;ht ('enter will manuallv 
map St'A fer eaeh satellite data tvpe He will use this information to evaluate the rela- 
tive aeeuraey of maiuial w'ver automated snow mapping teehmques for satellite data v>f 
different spatial resolutions 

I 4 SqiulffoAiuHf of thl« KeiNNin'h 

The immediate importanee would be in the improved utilitv of satellite datc\ for 
mapping S('A A redueed dependenee on more expensive hi^;h resolution data from the 
VSS or TV sei\sors would save oosts, and implementation of automated digital teoh 
niques would improve the effioienoy of state and federal water management prO|;rams 
In the loi\^ run. the impiwed understandii\^ of satellite data and their potentials and 
limitations, should provide broad benefits to both hydivlo^v and remote sensing; 


OF POOR QbAUT'i 
•7- 


Z ATMOSPHERIC RADIATIVE IHANSm MODEL ACX:0UNT1NG FOR BOTH THE ZENITH AND 
AaHUTH VARIATIONS Di THE RANATIVE FIEU) 

2.1. Synopns 

In remote sensing, the signature obtained by the satellite sensor is related to radi* 
ance, which is a function of elevation, location, and viewing angle, in addition to the 
astronomic and atmospheric conditions. 

In order to infer the radiometric properties of the medium underlying the atmo- 
sphere from remotely sensed data, atmospheric correction of the satellite data is 
necessary. In remote sensing for oceanography, quick atmospheric correction algo- 
rithms have been developed [Gbrdon, 1978; Gordon and dark, 1980] and successfully 
applied to the quantitative assessment of chlorophyll concentration in Southern Cali- 
fornia [5Vnif/i and Wilson, 1981], In these algorithms, the known spectral optical pro- 
perties of water are used to determine the aerosol path radiance at wavelength 670 
nr . where the subsurface upwelling radiance is negligible, atid then to estimate the 
aetosol path radiance at other wavelengths. However, in other wavelengths, or for 
other types of surfaces besides water, the application of these algorithms may be lim- 
ited, because of the different optical properties. Therefore, some general atmospheric 
correction algorithms are required. 

In our project, a different approach, which invokes an atmospheric radiative 
transfer model accounting for both the zenith and azimuth variation in the radiative 
field, is taken. The model deals with a plane parallel structured atmosphere, which is 
composed of different layers with each assumed to be homogeneous in composition and 
to have linear-in-r temperature profile. Moreover, the multiple scattering between 
such an atmosphere and the underlying medium (say, water, snow, or soil), which may 
have both the zenith and azimuth variation in radiative behavior, is taken into con- 
sideration. 

This approach is an extension of an azimuthally integrated atmospheric radiative 
transfer model developed by Wiscombe [1981], which is based on the results of earlier 
researchers [Qnmt and Hand, 1969] and his own investigations [Wiscombe, 1976a, 
1976b]. The additional parameter on azimuthal variation is handled by adding a new 
azimuth dimension into the computation space. In order to save computation time, 
the Fourier transform is performed over the azimuth domain. Then the resulting 
Fourier coefficients with different order stand for the original vector and work in a 
parallel manner, so that the array sizes for each order of coefficient in the computa- 
tion remain the same as for the azimuthally integrated case. The final radiance are 
obtained by Inverse Fourier transformation from the Fourier coefficients of those quan- 
tities, which are the results derived from the procedures similar to what is used for 
calculating the eizimuthally integrated radiance. 

The inputs of such a model are some astronomic and atmospheric parameters for 
each layer (or level), and the surface radiative property. Those parameters or proper- 
ties can be obtained either by some model or by measurements. The astronomical 
parameters are earth-sun distance and solar flux at the top of the atmosphere, and the 
atmospheric parameters include pressure temperature, cheniical composition of the 
air molecules, and the composition and size of the aerosol, water droplets, and ice cry- 
stals. The outputs of the model are the monochromatic radiemce and iiradiance of at 
each level. By integration, the total radiance and irradiance can be obtained. 

The potential application of this model in the atmospheric correction of remote 
sensing include several aspects. For example: after cedculating the upward -adiance at 
the top of the atmosphere for an atmosphere with black surface beneath ic, and then 
by removing this value from the radiance measured at the top of the atmosphere, one 
can get the modified radiance which accounts for the contribution of upward radiance 
at the surface level. 
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Also, this model can mimic the angular radiative field for a plane parallel layered 
atmosphere - earth system, in which the earth surface property is homogeneous from 
location to location. For low frequency spatial variation, this model may also offer good 
approximation. For more complicated cases with high frequency spatial variation, this 
model may approximate the average radiance over an area with certain size. The 
departure from the average may be used to derive more detailed information of the 
surface. 

Besides, this model will be useful to assess the other atmospheric correction algo- 
rithms, in which some aspects of the interaction between radiation and medium may 
be omitted. Because this model considers all the aspects of the interaction between 
photon and medium for a plane parallel structured absorbing and scattering system, if 
we apply the other algorithms on such a system, the departure of the results of the two 
models will give em indication of the performance of the algorithm to be assessed. 

2.2. Review of Radiation Model 

Once the atmosphere is subdivided into many plane parallel homogeneous layers 
according to some atmosphere model or some measurement data, the basic inherent 
radiative properties, such as the total optical depth, the single scattering albedo, and 
the phase function can be derived for each layer by Rayleigh scattering, Mie scattering, 
line absorption and continuum absorption theories. 

Based on those basic properties, the reflection and transmission functions, and 
the solar and thermal source function for each layer can be determined, and finally, 
the radiances at desired levels can also be cedculated, using the radiative transfer 
theory. 

In the following sections, we will deal with a A'-layer atmosphere plus underlying 
medium model. 

2.2. 1. RadiaHve Transfer Equation 

The radiative transfer equation for a plane parallel medium is given by [Chon- 
drasekhar, 1960] 

+ L{r,fx,y)) = J{T.fj,.<p). 

T is optical depth, and L{r.n,fp) is the radiance at level r along direction fi, ^ (where 
is cosine of zenith angle and ^ is azimuth angle). The source function is given by 

^ 2n 1 
n -1 


The last term represents an internal source. The internal source term in this case is 
similar to that for the azimuthally averaged case [ Wiscombe, 1976]. 

Qir.ti.ifi) = 

The thermal source is 

= (1-0)5(7’(t)) 

B{T{t)) is the Plank function at temperature T. where T is a function of optical depth. 
The direct source and specular "pseudo-source” arising from the diffuse - direct split- 
ting of the radiance are 


and 
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Q»ir.fi,<p) = (Mo) /’(t.m.PI-Ao.Po) 

Here (iq is the cosine of the solar zenith angle, Eo is the solar irradiance incident on the 
top of the atmosphere (normal to the beam), and p, is the directional specular 
reflectivity at the surface beneath the atmosphere. 


2.2.2. Interactixm Prineipla 

One way to attack this integro - differential equation is using the “interaction prin- 
ciple" [Ount and Ifunt, 1969] which, specified to a layer bounded by T/ and tj across 
wlUch 0 and phase function P are constant, assumes the discrete form 

L*(Tj) = t(Ti,Tj)L-{Tj) + t{Tj.Tj)L^Ti) + S"^(T/.Ty) 

^"(T/) = r[Tj.Tj)L*{Ti) + t{Ti,Tj)L-{Tj) + 2 -(T/.Ty) 


Radiance (t=T/,Ty) is a vector of mxn elements on a discrete angular space 

composed of m zenith and n azimuth angles 


(t) = 


L(t.±/x„^p,) 

L{f,±p.i,<P2) 




0<p,i< ■ ■ ■ <Hm^^ are a set of quadrature points for (0,1). 0s^i< • • • <^„<2 tt are 
equally spaced points in the interval 0-2rr. and t{Tj,Tj) are reflection and 

transmission matrices. E*(0,t) are internal source vectors. The r, t and S for each of 
the homogeneous single layer can be derived by some initialization scheme and the 
doubling method. 


2.2.3. InUialization 

In a recent letter, Wiscombe suggests that the Infinitesimal Generator Initializa- 
tion (IGI) is best [O^inf and Hunt, 1969]. It produces the t{Lt), f(Ar). and S*(At) for 
an optically very thin starting layer with optical depth At. After initialization, the 
r(T/.T/ 4 .i). f (t/.T/^i). and E*(t/,t/ 4 .i) for the entire homogeneous layer bounded by lev- 
els / and /+1 can be derived. 

The IGI claims 

t(At) = p*- c 

4tt 

f (At) = 7 - At p** C 

where I.M.C are matrices of size mn xmn. Af"* is the inverse of M. I is the identity 
matrix, M and C are diagonal matrices, 

M = [Pi 6ij <5«] 

C = [C* (5y dfct] 

and 

P** = [P{Pi.<Pk:^Pi.<Pk)] 

where liism. isjsm., liitin, and Kfsn, and pi's are the cosines of zenith angles 
derived by some quadrature rule (say, Gaussian quadrature rule), and the Ci's are the 
associated weights, and 6^ is delta function defined as = 1 if i =; , otherwise =0. 
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The source vectors for different source need to be handle separately. The IGl gives 
the general relation between source term Q‘s and the source vectors for the initial 
layer 


E*(At) = AT* 


/ 




Q{T,±tXi,<pz) 


dr 




The solar source vector for the initial thin layer is 

« («T) = ^ (' - 

and the thermal source term 


However, the thermal source vectors for the entire homogeneous layer with linear-in- 
tau temperature prohle within it can be derived in a quite simple fashion without 
directly using S^(At). 


2.2.4. Adding and Doubling Method 

By applying the interaction principle to two adjacent layers, the reflection and 
transmission matrices and the source vectors for the combined layer can be derived if 
the corresponding quantities are known for each of these two layers [Oanf and Hunt, 
1969). 

Consider two adjacent layers bounded by planes x.y.z. By the interaction princi- 
ple, we have 

i*(i) + r(z,y) L-{y) + V{x.y) 

L~[x) = r{x.y) L*{x) + t{x.y) L'{y) + l.'{x.y) 

L*{z) = t{z,y) L*{y) + r(y.z) /.‘(z) + Z*{y.z) 

L~{y) = r{y.z) L*{y) + t{y.z) L'^z) + E'(i/.z) 

Since z.y.z are entirely arbitrary, we can use the interaction principle and v/rite 

L*(z) = t{z ,x) L^{x) + r{x,z) L~{z) + E*(z,z) 

L~{x) = t{x,z) L*{x) + t{x.z) L~{z) + E"(z,z) 

These equations can be obtained by eliminating L^{y) and L~{y) from those before. 
Then the following relations hold 

r(z,z) = r{y,x) + t{x.y) [/ - r{z,y) r(z,i/)]'‘ r{z,y) t{y.x) 

r(z,z) = r(y,z) + f(z.v) [/ - r{x ,y) t{z ,y)Y^ r{x,y)Hy,z) 

t{z.x) = t{z.y) [I -r(z,v)r(z,v)]-' f(y.z) 

f(z,z) = f(z,y) [/ -r(z,y)r(z,y)]-‘ f(y,z) 

E* (z,z) = f(z,y) [/ - r(z,y) r(z.y)]'* r(z.y) E"(y,z) + 

f (*.y) [/ -T-(*.y)’'(*.y)]'‘ s*(y,z) + E*(z,y) 

E'(z.z) = f(z,y) [/ -r(z,y) r(z.y)]-' r(z.y) E^(y.z) + 

t(x,y) [/ -r(z,y) r(z,y)]'* E“(y,z) + E'(z.y) 

These are the formulae for the so called "adding" method. If the two layers are 
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identical, the method is called "doubling." If the Initial layer is chosen such that 
At = (t/», - t/)/2^ 

where TV is a integer, and (t/ 4 ,| -t/) is the optical depth of the 7+1 th layer in the 
multi-layer system, then the reflection and transmission matrices and source vectors 
for such a homogeneous thick layer can be built up quickly by "doubling" TV times; 

= ffi (7 -T„ r„)"‘ 

Tfi^l ~ Tn + fn (7 ~ T|j Tu) * 
r„=r(2"AT), and OaJniAT 

Since neither of the internal sources, including the solar and thermal sources, are 
consteint with optical depth, nor the changes of them with optical depth the same, they 
need to be treated separately. For the solar source, the doubling scheme takes a spe- 
cial form: 

~ ffi Tf> (Cu Tjj E^.n + Ed.n) ^ *n 
^e.n + l “ (Tfi *« Sd.n) ^ 

For initial layer, n=0, eo=e”*^^^, and e„ 4 .,=e^. 

For a layer with linear-in-r temperature profile, the doubling method gives the 
final resulting thermal source vectors as [ fHscom.be, 1976]; 

+ Bn) Yn ± B’ Zn 

where 

^o = ^[nT/)] 

Bn = B[T{t,,,)] 

B' = [(£^JV-^o)/(T/tl-T/) 

Zs can be obtained recursively: 

^n*\ “ 2'n + Pn Y^ + r„) (Zf^ — Y^) 

Pn = 2 "-*At 

Yn,, =(tn r„ r„ + r„ + 7) 

with the initial values Zg=0, Po=)JAt and Fo=(1~®)At M~^U . U is a. vector of ones, and 
r=(7-r„r„)-« 

Then the total source vectors for the entire homogeneous layer are the sum of all 
the different types of internal sources. 

= ^i.N + + • • 


2.2.5. Calcxdation of the Internal Radiance 

Having known the reflection and transmission matrices, and the source vectors for 
each layer in the multi-layer system, we can build the internal radiance field in the 
atmosphere by the adding method. Using formulae of the Interaction principle, we 
have a set of simultaneous equations: 

7-*(T/*|) = f (T/ + ,.T/) L*(T/) + r(T/.T;*,) + + S*(T/,T/h) 

7.'(t/) = t(T/*,.T/) L*(ti) + f (T/.T/+,) L'iTj^i) + 

for fts7<7C. where K is the total number of layers in the system, and To= 0 and is 
total optical depth of the atmosphere. In addition, two sets of boundary conditions 
need to be saUjfled: 
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L-(tk) = /fc L^(tx) ^ c ff(Tc) * 


is the bottom boundary condition. Jfc i> the surface diffuse reflection matrix, fr{fiQ) is 
the BRDF (bidirectional reflectance-distribution functions) vector for the incident 
beam, and To is the temperature of the surface. 

It is convenient to solve this set of equations in the following recursive manner. 

1. Compute recursively the intermediate partial interactive quantities from top to 
bottom for each layer (or level). 

This is so called forward pass. In this stage, the calculation starts at the top layer, 
then goes down, and for each layer only the contribution made by the internal source 
of the current layer and by the radiation coming from above, and the interaction 
between the current layer and the layers above are taken into account. There are six 
quantities calculated. They represent the partial upward and downward radiance 
fields, and incident-from beneath reflection and transmission operator, respectively. 
The precise physical meanings of these quantities can be understood from the following 
recursive formulae and the initial values for them. 


Vf = [/ - r(T;,,.T/) /?£)-,]-» [I:-'(t/,T/*,) + r(T/*,T;) Vf., ] 

V,* = Vf., + RE,., V,- 

K/= £^(T/,T/+i) + f(T/+,,T/) V,* 

H, = [/-1-(T/*,,T/) f(T/.T/+,) 

G, = RE,., H, 

RE, = r(To,T;+l) = r(T/,T/t,) + f (T/h.T/) G, 
where 7=0, • • • ,K-\ and the initial values for them are given by 
Vo" = £'(to,t,) + r(T,,To) L*(to) 

Vi = L*(to) 

Vl = I:*(to,t,) + f (tj.to) L*(tp) 

RE„ = r(To,T,) 


Co = 0 

77o = t (to.Ti) 


2 Calculate the radiance at the atmosphere - earth interface. 

From the forward pass we have calculated the total contribution of the internal 
sources and the radiation coming from the top boundary to the downward radiance at 
the surface level, which is Vj,., . Using the interaction principle again, we can come up 
with the final result for the downward radiance at the earth surface. 

L*{tk) = [I -REK.,RcY^[rK., + 

REk., ( + c B{Tc))] 

The upward radiance at this level can be obtained from the lower b'^undary condition. 

3 Compute recursively the downward and upward radiances from bottom to top. 
This is so called backward pass. 
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. successively for the combined layer, each time adding one more layer. Finally, obtain 

the corresponding quantities for the entire atmosphere. 

Then, calculate the upward radianc; using 

^■(to) = k(Tjr.To) + f (to.tjt) Rq {I RgY'^ f (T/T.To)] L*{Tq) + 

S^'(to.Ta-) + t (to.t^) Rc (/ - r (tq.Tjt) rc)-> x 

[S^*(To.Tjr) + r(To.Tjr) [tB{Tc) + /r(/io)]] 

L'.3. Coaxmenta 

There are several advantages with this model. 

1. It considers all aspects of the radiation - medium interaction, including multi- 
ple scattering occurring within the atmosphere and between the atmosphere and the 
underlying medium. The crymplete consideration of ail the aspects in multiple scatter- 
ing is especially important when the surface has a high reflectance. 

2. For surfaces with moderate area size, the algorithm can be used to model the 
radiance field above it. The modeling accuracy will be good as long as the plane paral- 
lel configuration is valid for the situation. 

I 3. This algorithm can also be used to calculate the averaged atmospheric correc- 

tion amount. When this amount is removed from the remotely sensed data, the result- 
ing radiance values will be good indication of the local surface variation. 

4. The algorithm will also be usef'^ for removing cloud interference, as long as the 
clouds are not opaque in the wav ilength concerned. 

I 5. Multiple sources, including thermal source, direct solar beam source, and pos- 

I sibly even specular ^seudo source, are taken into account. 

[ 6. If the upward radiances at the top of the atmosphere are measured systemati- 

cally at a set of viewing angles, then the angular radiation field at the surface level can 
be inferred by inversion. 

i The disadvantages of this algorithm include the relatively high computation time 

consumption and relatively large number of input parameters. 

3. SOLVING THE ATIiOSPHERIC CORRFXniON PROBl£U FOK AN INHOMOGENEOUS SUR- 
FACE USING A MONTE CARLO METHOD 

3.1. IntroductioD 

In remote sensing, Landsat satellites look at the eeu*th at near nadir position and 
with a very narrow range of viewing angles. Therefore, by fair approximation for our 
present atmospheric correction algorithm. Landsat looks straight downward at tlie 
earth surface. 

According to the interaction principle [Orant and Hunt, 1969]. the upwelling radi- 
ance for a layered structure atmosphere at the height of the satellite can be 
expressed. 

L"(to) = r(T,.To) L*{tq) + f (to.t,) L‘(t,) + S"(to,t,) 

To=0 if the optical depth at the satellite location, and Tg is the optical depth at the 
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ground level. The layer concerned is bounded by the levels at tq. and t,. r(T,.To) is 
the reflection matrix of the atmosphere of the layer of atmosphere when the incident 
light L*{to] comes from the top of the layer, f (tq.Tj) is the transmission matrix when 
the incident light Z "(tj) comes from the bottom of the layer, and £"(to.t,) is the con- 
tribution of the source within the layer to the upwelling radiance at the satellite level 
To 

Note that Z*(to), L~{rg). and £"(to.t,) are vectors describing the angular distribu- 
tion of the radiauice, while t(t,.To) and t{To,Tg) are square matrices of corresponding 
size. Also, note that the incident radiances at boundaries are independent of location 
by the above formula. 

For a given atmosphere, the terms r(Tp.To) L*{tq) and S“(to.Tj) can be calculated 
from an atmospheric radiative transfer model accounting for both zenith and azimuth 
dependence in a absorbing and scattering layered atmosphere. 

The element of Z’(to) at the nadir viewing angle, denoted as L“(To.d„.9Pv). can be 
obtained by the satellite measurement. The main concern in remote sensing is to get 
the quantity L"(Tg.d„.^v). which is an element of Z“(Tj) and where and (py are the 
zenith eind azimuth angles of the viewing direction of the satellit j. For an optically thin 
atmosphere, the term is the leading one contributing to the signature 

Z'(to.'5v.^v)- It the pattern of the surface radiance is known or a Lambertian property 
is assumed, the surface upwelling radiance can be derived given the satellite measure- 
ment and the atmospheric profile, using a layered atmospheric radiative transfer 
model. 

However, for most remote sensing, the situation is more complicated, since the 
upwelling radiance at the earth surface is not independent of location. We can still 
consider the atmosphere as a layered medium, and within each layer homogeneity can 
be assumed. In this case, the terms corresponding to r(Tj.To) Z*(to) and E"(to.t,) £ire 
still independent of location and remain the same as for a homogeneous surface. 
Theiefore, these two terms can be calculated from the layered radiative transfer 
model assuming the underlying surface is black and non-emitting. The only problem is 
that the contribuHon of the surface upwelling radiance to the upwelling radiance 
detected by the satellite sensor cannot be simply obtained from the term 
f (to.Tj) because L~[jg) is now location dependent and the corresponding 

transmission function is also different from the term 1 (to,Tj) which is related to the 
integration of the contribution of ground to the radiance at the satellite level over an 
infinite homogeneous surface area 

Because we are dealing vrith 2-dimensional surface, we can express the satellite 
measurement by the following equation; 

Z'(To,1>v.5<’v.f.’7) = Arfm(To.t>v.<<’v) + AZ--(Ti.iJv 

yu = cosiJ„ and z ,y are the coordinates for the ground position for the measured pixel, 
euid i,T] are the horizontal coordinates at the satellite level. The first term at the right 
hand side is the element of AJin(To) at viewing angle and the latter is obtained 

from the layered structure atmospheric radiative transfer model. 

Zaim(To) = t(t,,To) Z*(To) + E"(To.T,) 

If we can estimate then we can obtain Z~(t,.i>v.^v^^) using 

remote sensing measurement and the layered atmospheric model. The Monte Carlo 
method is one of the suitable candidates for this purpose. 
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3.2. Monte Carlo Apiroaeh 

Only the transmission problem is treated, because the reflection and source prob- 
lem can be solved by the previous model. For a system composed of a layered atmo- 
sphere and an underlying inhomogeneous surface, a combination of the azimuth depen- 
dent radiative transfer model for the term and the Monte Carlo method for 
transmission of surface upwelling radiance throngh the atmosphere layer can give a 
complete expression of the satellite measured radiance. 

The proposed model is based on the fact that the contribution of the surface 
upwelling radiemce to the signature of a certain pixel is made by not only the ground 
element corresponding to the pixel concerned, but also the surrounding area as well. 
By the reciprocity principle, the pattern of the contribution made by the ground area 
can be mimicked by a reverse process, a process in which a beam of photons impinges 
at a given point at the top of the layer and flneilly some of them get to the surface and 
make a spread pattern, the so called "point -spread" function, in other words, if we can 
find the point-spread functi 5n for a beam impinging on a given atmosphere as a func- 
tion of incident angle, we can obtain the quantitative expression for the pattern of the 
contribution of surface upwelling radiance to the satellite measurement. This relation 
can be written as the following formula: 
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The IS the basis of the present algorithm, where t is 




/.(T2.X2,y2,'dg.^2) cosiJa duz 


AL is the increment of radiance at location (T,.x,.y,) and in the direction ("di.^Pi) con- 
tributed by the radiant flux at location (Tg.Xj.yg) in the direction 

Note that for this process, the conservation of photons docs not hold because 
some of the incident photons may be absorbed and some of them may get out of the 
layer through the upper boundary. In either case, the photons are lost and have noth- 
ing to do with the transmissi poirt-spread function. Besides, in this configuration, we 
need not consider the emissit .* problem because the previous adding model has taken 
care of that process already, euid this makes our Monte Carlo method relatively simple. 

We need to determine the point-spread function using the Monte Carlo method. 
Several conditions are given: 

1 . The beam of photon impinges at the top of the layer of atmosphere concerned. 

2 The atmosphere has layered structure and the profiles, including molecular 
absorption and scattering. Mie scattering components, and temperature, arc known. 

3 f’or each sublayer, homogeneity is assumed with the scattering and absorption 
properties assumed to be the average value derived from the given profile. 

4 The relation between optical depth and height should be derived from the 
atmospheric profile, and be used throughout this algorithm, since this knowledge is 
crucial for calculating the horizontal position at which the photon gets out of the atmo- 
sphere 

The essence of the Monte Carlo method is that the scattering and absorption of the 
photon can be statistically simulated by a sequence of photons randomly colliding with 
the medium particles before they finally get absorbed or escape from the medium con- 
cerned Hach scattering or absorption is a random event of collision. However, the 
general trend is governed by some statistical rule according to some probability func- 
tion of the processes. 

In order to mimic this random process, we need some quantities based on some 
statistical rule. These are the distance traveled by the photon between two random 
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collision events, the direction of each traveling path, and the chance for the photon to 
be absorbed upon collision. Also, we need to keep track of position of the photon in 
three-dimensional space in order to know whether the photon is in the atmosphere, in 
which sublayer the photon is, and how far it is from the original horizontal position. 

The following sections give the mathematical expressions of these events and 
quantities. Similar descriptions c«ui be found in some representative papers [ Cash-well 
and Fhjerett, 19&9; House and Avery. 1969; Pearce, 1977], 

3.2.1. PYee f\ith Length 

The distance traveled by a photon between collisions is called the free path length. 
When It is measured in the seune manner as the optical depth, it is dimensionless and is 
called optical distamce The probability density function of the occurrence of noncolli- 
sion for an optical distance I is 

P(0 = e ‘ 

The probability that no collision occurs in the range of optical distance from 0 to f is 
1 

- J p{l ) dl' = \ - e~^ 

0 

This equation sets up a unique relation between a random number r, (0.CKr,<1.0) and 
an optical distance I , and I can be solved explicitly as 

I = -!n(l -r,) = |ln(l -r,)| 

Therefore, the distance between scatterings for each path of each photon can be deter- 
mined according to a random number r, generated in some manner. 

3.2.2. Chance of Absorpt-ion 

The single scattering albedo 0 of the current sublayer gives the probability 1-S) 
that the photon is absorbed upon collision. Conceptuedly, we can use the idea of photon 
bundle instead of single photon. Suppose initially a photon bundle impinges into the 
atmospheric slab. After collision, only a fraction of the photon bundle keeps going and 
a portion is absorbed F’or each collision, the current fraction is multiplied with the 
single scattering albedo of the current sublayer to get the fraction of the photon bun- 
dle continuing Such a process is repeated until the remaining portion of the photon 
gets out of the slab or is less than some threshold. 

3.2.3. Direction of Scattering 

F’or each scattering event, two independent angular variables can be obtained 
from random processes, the scattering angle 0 and the azimuth angle which is meas- 
ured in the plane perpendicular to the original travel direction 

The scattering angle 0 is determined in the following way: 
define 

e 

P(0) = /p(0')sin0'd0' 

0 

0' is a dummy veiriable, p(0’) is the phase function for scattering angle 0' for the 
current scattering sublayer. Because the integration of p(0) over the range from 0 to 
TV is 2. we need to multiply the quantity expressed by a factor of O.b to normalize it. 
Then we can relate such a quantity to a random number rg generated for the random 
scattering process to determine the scattering angle 0. 

^P(0) =ra 
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For 0fi'»J<7T/ 2, 

fz~f\ (<=“■*> 

For iJ=n/ 2, 

/2 = 0 

T* is the total optical depth. The horizontal position of escaping in the manner 
described earlier. The remaining portion of the photon bundle fi~fz will still travel 
within the slab. The distance traveled related to random number r, can be determined 
from 

I = |ln(l-r.(l-=^)| 

/ 1 

Such a process along with the absorption process upon collision for non-conservative 
scattering makes the remaining portion of the photon bundle smaller amd smaller amd 
eventually the remaining fraction is so small that it can be neglected. Such a treat- 
ment can save computation time for given precision. 

3.2.7. Piscretization of Armies and Distance 

When the photon travels within the slab, the locations of the photon are calculated 
in a continuous manner so that no error is accumulated. However, at the point the 
photon escapes at the bottom of the slab, we have to record the escaping location and 
direction of the photon in a discrete manner. Four dimensions of discretization may 
be applied since the distribution of the escaping photon is a function of location coordi- 
nates of I ,y and angular coordinates An £U*ray of 4 dimensions takes a lot of com- 
puter memory space, and to make the recorded photon number in each element of 
such an array statistically meaningful, the total number of photons needed in the 
experiment must be tremendously large. Therefore such an algorithm is infeasible. 
Under the assumption that the photon impinges at the top of the slab perpendicularly, 
the point-spread function is axisymmetric and we need to record the information along 
a single radius only instead of over the entire x.y plane. Therefore a recording array 
of only 3 dimensions is needed. We fix the single recording radius along the x axis and 
put on it the results from the location away from it by rotation. The rotation angle (Pr 
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tan(fl, = yz/xz 

Note that the signs of Xz and yz allow (pr to be determined over the range 0 to 2rr. After 
rotation, the zenith angle remains the same while the azimuth angle is changed. 
Therefore the new zenith angle iJ is 

= T?2 

and the new azimuth angle <p is 

tp =<pz~pj. 

The distance of the total horizontal displacement is 
r = + yz^ 

Later on we will see that for a Lambertian surface, the discretization of and p are not 
necessary, only the spatial discretization is needed. For such a case, the computer 
resources needed can be reduced, and the non-axisymmctric case would be computa- 
tionally possible. 
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For the axlsymmetric case, the discretizations of r, ‘d, and (p can be performed in 
the following m&aner: 

For the r -dimension, an equal increment discretization scheme is carried out. 
The increment of r should be equal to or less than the ground resolution of the pixel. 
For zenith angle t>-dimension. the double Gaussian quadrature rule for cosine of zenith 
angle is the suitable one because the result obtained will be comparable to those pro- 
duced in the layer adding radiative transfer model. For the azimuth angle domain, the 
best discretization is the equal increment scheme. The range covered by each element 
can be determined by finding the midpoints between the adjacent discrete points. 


3l2.B. Calculation of the Point Spread Function 

In digital image processmg, we need to express the point-spread function {PSF) in 
matrix format. If we superimpose such a matrix over a matrix whose elements 
represent surface irradiance. then the product of two corresponding elements in these 
two matrices represents the contribution made by the irradiance from that particular 
ground element to the radiance measured by the satellite, and the sum of the products 
within the PSF window is the total contribution made by ground radiation to the signa- 
ture of the pixel that coincides with the central element of the PSF matrix. 

In order to construct such a matrix, we need some corrections to the results 
obtained directly from the Monte Carlo method. 

First is the correction for the difference between areas covered for a fixed incre- 
ment of radius Ar but for different values of radius r. This is obvious, since the central 
point’s influence ranges from 0 to and therefore the influence area is (n/4) Ar^. 
Other points’ influences range from r-)5Ar to r+JJAr with area 2rrr Ar. Therefore it is 
necessary to make the area correction to get the contribution per unit area. 

Secondly, the difference between the area of a circle and the area of square or 
rectangle must be corrected. For a matrix, the most convenient shape for each ele- 
ment is square or rectangle. 

After these corrections, the contribution made by the surface radiation from a 
unit area of Lambertian surface for the central ground element emd the surrounding 
elements can be expressed as 


Wq 


S i: C 

¥ ^ 


Nr 


ti;, = 


S D N[r.A<p) C 

f V 



u>a is the central element of the matrix of the point-spread function, and -uv is an ele- 
ment other than the central one with distance fror- center r. N is the 3-dimensional 
array obtained in the Monte C£U"lo algorithm, and Nj m the total number of photons 
incident at the top of the atmospheric slab in the experiment. C corrects for the 
difference between the area of a circle and that of a square or rectangle. For a square. 
C = 4/ TT. From these equations, it is obvious that only the radius dimension is neces- 
sary for Lambertian surface. 

The convolution of such a point-spread function matrix with the matrix of ground 
radiance, in which each ground element is Lambertian, gives the total contribution of 
surface upwelling radiation field to the satellite measurement, and the satellite meas- 
ured radiance can be expressed as 

L (T(j,iJv >^ j) ~ Ao<m(To, »?v 
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h is the point-spread function matrix of size it is convolved with a ground matrix 

of size NxN. The elements of h outside the MxM window size are all zero. The last 
term on the right hand side is the convolution term. 

The point-spread function matrix thus defined only works for the Isotropic surface. 
For anisotropic surface, we need to define the contribution factors for the central 
ground element and the surrounding elements, respectively. 


A/(T,.xo.yo) Nt 


= :r^E E ETA{Tg.i^.<p.Xo,yo) N(ro.i>.5fl) 
E E C At 

A.. - J-1 

^ Q M{Tg.x^.y^)r^j Nt 


ETA is the anisotropic factor for upwelling radiation under given illumination condition, 
and is defined as 


ETA{Tg.^,^.x,y) 


nL (Tj.d.^g.i.y) 

M{Tg,x.y) 


M is the exitance at given location. Then the satellite measurement of the radiance is 


+ £ hJ^{x^^.y^.j) M{Tg,x^.yj) 

<s0 >s0 

IS the element of the matrix composed of the contribution /ly’s. For a Lambertian 
surface, a simpler formula can be used. 

* t t MX. -i.yi-i) L~iTg.Xi.yj) 

t=0 j=0 

h is the element of the matrix composed of weight ioq and ti^'s. 

From these equations, it is obvious that if the point-spread function is derived 
from the Monte Carlo method and if the surface exitance and the pattern of the radi- 
ance angular distribution at each surface location are known, the vertically upwelling 
radiance at satellite level can be determined. 


This is the forward aspect of our problem. However, in remote sensing, the 
reverse aspect of the problem is of more significance, but it is more difficult, too. In 
the following section, we discuss the inverse problem. 
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3.3. Retrieval of the Surface Exitance 

After removing the atmospheric reflection and source term from the satellite 
measurement, we have an image g. Also, we can construct a matrix H from the point- 
spread function such that the size of H is comparable to the lexicographic form, by 
which the rows of a matrix are stacked to form a vector of the matrix y. and the lexico- 
graphic form of the corresponding matrix of surface exitance / . The relationship 
between / , y , and H is 

g=H f 

If there is no noise other than the atmospheric effect, then the surface exitance can be 
retrieved by the inverse filter given y . 

/ = ^-‘ y 

The Fourier tremsform technique can reduce the computation time for the inverse 
filtering [i4ndretws and Hunt, 1977]. In the presence of noise, other types of filters may 
be employed using the point-spread function along with some a priori knowledge of 
noise [Andrews and Hurd, 1977]. However, those procedures are more complicated 
and the extent to which the result can be improved strongly depends on the 
signal /noise ratic. The Fourier transform technique similar to that used in inverse 
filtering can also used in those filters, to reduce computation time. 
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